 rm(list = ls())

sel=list.files(pattern=".snp")


sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49")
sel1=paste(sel1,".snp",sep="")
countdiff=c()
for (i in 1:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="T",]
if(dim(file)[1]!=0){
file$tumor_var_freq=as.numeric(gsub("%","",file$tumor_var_freq))/100
file$normal_var_freq=as.numeric(gsub("%","",file$normal_var_freq))/100
file$total_reads=as.numeric(file$normal_reads1)+as.numeric(file$normal_reads2)+as.numeric(file$tumor_reads1)+as.numeric(file$tumor_reads2)
sels=which(file$total_reads>=10)
file=file[sels,]
file$DU=abs(file$normal_var_freq-0.5)
file$DA=abs(file$tumor_var_freq-0.5)
file$p.diff=lapply(1:nrow(file),function(x) fisher.test(rbind(c(file[x,5],file[x,6]),c(file[x,9],file[x,10])))$p.value)
file$p.biasA=lapply(1:nrow(file),  function(x) binom.test(c(file[x,9],file[x,10]),p = 0.5)$p.value)
file$p.biasU=lapply(1:nrow(file),  function(x) binom.test(c(file[x,5],file[x,6]),p = 0.5)$p.value)
bias=(file$DU>=0.25 & file$p.biasU<0.05) | (file$DA>=0.25 & file$p.biasA<0.05)
diff=(file$p.diff<0.05)
count=paste(file[,1],file[,2],file[,3],file[,4],sep=":")
countdiff=unique(c(countdiff,count[which(bias & diff)]))
}
}

rt=data.frame(as.matrix(countdiff,1,length(countdiff)))
names(rt)=c("unitID")
rt=tidyr::separate(rt,unitID,into=c("chr","pos","ref","var"))
rt$unitID=paste(rt$chr,rt$pos,sep=":")
positions=as.character(unlist(rt$unitID))

i=1
file=read.table(sel[i],head=T,sep="\t")
file=file[file$exits=="T",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],":",file[,2],sep="")
rownames(file)=file$unitID
final=file[positions,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final=final[complete.cases(final$unitID),]
final$tumor_var_freq=as.numeric(gsub("%","",final$tumor_var_freq))/100
final$normal_var_freq=as.numeric(gsub("%","",final$normal_var_freq))/100
final$total_reads=as.numeric(final$normal_reads1)+as.numeric(final$normal_reads2)+as.numeric(final$tumor_reads1)+as.numeric(final$tumor_reads2)
sels=which(final$total_reads>=10)
final=final[sels,]
final=final[,-8]
final$DU=abs(final$normal_var_freq-0.5)
final$DA=abs(final$tumor_var_freq-0.5)
final$p.diff=lapply(1:nrow(final),function(x) fisher.test(rbind(c(final[x,1],final[x,2]),c(final[x,4],final[x,5])))$p.value)
final$p.biasA=lapply(1:nrow(final),  function(x) binom.test(c(final[x,4],final[x,5]),p = 0.5)$p.value)
final$p.biasU=lapply(1:nrow(final),  function(x) binom.test(c(final[x,1],final[x,2]),p = 0.5)$p.value)

colnames(final)=c(paste(strsplit(gsub(".snp","",sel[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp","",sel[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID",paste(gsub(".snp","",sel[i]),c("DU","DA"),sep="_"),paste(gsub(".snp","",sel[i]),c("p.diff","p.biasA","p.biasU"),sep="_"))
}

i=3
for (i in 2:length(sel)){
file=read.table(sel[i],head=T,sep="\t")
file=file[file$exits=="T",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],":",file[,2],sep="")
rownames(file)=file$unitID
final1=file[positions,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final1=final1[complete.cases(final1$unitID),]
final1$tumor_var_freq=as.numeric(gsub("%","",final1$tumor_var_freq))/100
final1$normal_var_freq=as.numeric(gsub("%","",final1$normal_var_freq))/100
final1$total_reads=as.numeric(final1$normal_reads1)+as.numeric(final1$normal_reads2)+as.numeric(final1$tumor_reads1)+as.numeric(final1$tumor_reads2)
sels=which(final1$total_reads>=10)
final1=final1[sels,]
final1=final1[,-8]
final1$DU=abs(final1$normal_var_freq-0.5)
final1$DA=abs(final1$tumor_var_freq-0.5)
final1$p.diff=lapply(1:nrow(final1),function(x) fisher.test(rbind(c(final1[x,1],final1[x,2]),c(final1[x,4],final1[x,5])))$p.value)#计算各种统计值
final1$p.biasA=lapply(1:nrow(final1),  function(x) binom.test(c(final1[x,4],final1[x,5]),p = 0.5)$p.value)
final1$p.biasU=lapply(1:nrow(final1),  function(x) binom.test(c(final1[x,1],final1[x,2]),p = 0.5)$p.value)
colnames(final1)=c(paste(strsplit(gsub(".snp","",sel[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp","",sel[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID",paste(gsub(".snp","",sel[i]),c("DU","DA"),sep="_"),paste(gsub(".snp","",sel[i]),c("p.diff","p.biasA","p.biasU"),sep="_"))
final=merge(final,final1,by="unitID",all=T)

}
}
final2=data.frame(lapply(final,as.character),stringsAsFactors=FALSE)